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We report a new implementation for axisymmetric simulation in full general relativity. In this 
implementation, the Einstein equations are solved using the Nakamura-Shibata formulation with 
the so-called cartoon method to impose an axisymmetric boundary condition, and the general rel- 
ativistic hydrodynamic equations are solved using a high-resolution shock-capturing scheme based 
on an approximate Riemann solver. As tests, we performed the following simulations: (i) long-term 
evolution of non-rotating and rapidly rotating neutron stars, (ii) long-term evolution of neutron 
stars of a high-amplitude damping oscillation accompanied with shock formation, (iii) collapse of 
unstable neutron stars to black holes, and (iv) stellar collapses to neutron stars. The tests (i)-(iii) 
were carried out with the F-law equation of state, and the test (iv) with a more realistic parametric 
equation of state for high-density matter. We found that this new implementation works very well: 
It is possible to perform the simulations for stable neutron stars for more than 10 dynamical time 
scales, to capture strong shocks formed at stellar core collapses, and to accurately compute the mass 
of black holes formed after the collapse and subsequent accretion. In conclusion, this implemen- 
tation is robust enough to apply to astrophysical problems such as stellar core collapse of massive 
stars to a neutron star and black hole, phase transition of a neutron star to a high-density star, and 
accretion-induced collapse of a neutron star to a black hole. The result for the first simulation of 
stellar core collapse to a neutron star started from a realistic initial condition is also presented. 
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I. INTRODUCTION 

In the 1980s, one of the most important issues in the field of numerical relativity involved performing simulations 
of rotating stellar collapse with the assumption of axial symmetry. Simulations of rotating stellar collapse in full 
general relativity were first performed by Nakamura and collaborators [1,2]. Using the (2-|-l)-|-l formalism developed 
by Maeda et al. [3], they succeeded in performing simulations of a rotating collapse of massive stars to black holes. 
They used cylindrical coordinates [w, z) with a nonuniform grid spacing and with at most (42,42) grid resolution for 
{vj, z) because of restricted computational resources. 

To compute gravitational waves emitted during gravitational collapse to black holes, Stark and Piran [4] subse- 
quently performed simulations similar to those of Nakamura et al. , adopting spherical polar coordinates with a typical 
grid size (100, 16) for (r,9). The distinguishing feature of their work is that they adopted the Bardeen-Piran formal- 
ism [5], which is well suited for computation of gravitational waves in the wave zone. As a result of this choice of 
formalism, they succeeded in computing gravitational wave forms, and clarified that the wave forms are characterized 
by the quasinormal mode of rotating black holes formed after gravitational collapse and that the total radiated energy 
of gravitational waves is at most 0.1 % of the gravitational mass of the system [4]. 

Since the completion of their works, no new work in this field was done for the next 15 years [6]. Although several 
questions that they originally wished to answer have been answered by their simulations, it was not feasible to perform 
sophisticated astrophysical simulations in the 1980s. This is likely due to the fact that the computational resources were 
severely restricted, and, in addition, that techniques in numerical relativity such as methods to provide realistic initial 
conditions and to perform the long-term simulations were not sufficiently developed. As a result, there still remain 
many unsolved issues in astrophysics and general relativity that can be studied with axisymmetric hydrodynamic 
simulations in full general relativity. Among them, realistic simulations of rotating core collapse of massive stars, 
which thereby become black holes or protoneutron stars in full general relativity, have not yet been performed. Stellar 
collapse is a common phenomenon in the universe, and hence, understanding the formation mechanism of black holes 
and neutron stars in nature is one of the most important issues in astrophysics. Actually, the study of the formation 
of rapidly rotating black holes with surrounding accretion disks in stellar core collapse is currently one of the hot 
topics in connection with a hypothetical scenario for the central engine of 7-ray bursts [7]. To date, simulations of 
a rotating collapse of a massive stellar core have been done in the Newtonian gravity [8-13] or in an approximate 
general relativistic gravity [14] using the so-called conformal flatness approximation (or the Isenberg- Wilson-Mathews 
approximation). In rotating stellar core collapses, general relativity plays an important role. As demonstrated in [14], 
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general rclativistic effects modify the collapse, bounce, and amplitude of gravitational waves emitted significantly, 
even in the formation of neutron stars. Of course, general relativity plays a crucial role in the formation of black 
holes. Thus, general rclativistic simulation is inevitable to precisely understand the nature of stellar core collapses. 

One long-standing issue for axisymmetric simulations in full general relativity has been to develop methods in 
which the accuracy and stability for a long-term simulation can be preserved. In axisymmetric simulations, we have 
in general used cylindrical and/or spherical polar coordinate systems, which have coordinate singularities at the 
origin and along the symmetric axis raj = 0. At such coordinate singularities, the finite differencing scheme has to 
be changed, resulting often in numerical instabilities. To stabilize computation, we have often been required to add 
artificial viscosities around the coordinate singularities to stabilize the numerical system [15]. 

Recently, the Potsdam numerical relativity group has proposed the so-called cartoon method by which a robust 
numerical relativity implementation for axisymmetric systems can be made [16]. The essence of their idea is that 
the Cartesian coordinates {x, y, z) could be used even for simulations of axisymmetric systems if the Einstein field 
equations are solved only for the y = (or x = 0) plane, using the boundary condition at y — ±At/ (or x ~ ±Aa;) 
provided by the axial symmetry. (Here, Aa;, Ay, and Az denote the grid spacing.) Since the field equations are written 
in the Cartesian coordinate system, we neither have singular terms nor do we have to change the finite differencing 
scheme anywhere, except at the outer boundaries. Thus, it is possible to perform a stable and accurate long-term 
simulation without any prescription or artificial viscosities, but only by a minor modification of a three-dimensional 
implementation that has already been developed [17-20]. 

Other important progress has been made regarding computational resources. Current large-scale supercomputers 
that we can use are typically of several hundred Gbytes memory. Necessary memory in an axisymmetric simulation 
with double precision, with N'^ grid points, and with Ny variables is 

— '*<^)'(^). 

where Ny is ^ 200 in our general rclativistic implementation. This implies that the memory of current supercomputers 
is large enough to carry out an axisymmetric numerical simulation with N ~ several thousands. Using N of order 
10'^, it is feasible to carry out a wcU-rcsolved simulation and a careful convergence test, changing the grid resolution 
for a wide range from N ~ 100 to 1000. This situation is in contrast with that of three-dimensional numerical 
relativity, since it is still very difficult to carry out a three-dimensional simulation with N ~ 10^, for which the 
required computational memory is of order TByte. 

Motivated by the status mentioned above, we recently started a project in axisymmetric numerical relativity. In 
[21], we reported a numerical hydrodynamic implementation for axisymmetric spacetimes that is made using the 
cartoon method and incorporating a hydrodynamic implementation in the cylindrical coordinates. In that paper, we 
presented numerical results for simulations of rotating stellar collapse adopting simple initial conditions and simple 
equations of state to investigate the effects of rotation on the criteria for prompt collapse to black holes in an idealized 
setting. We demonstrated that the axisymmetric implementation and current computational resources are well suited 
to systematically perform stable and well-resolved hydrodynamic simulations in axisymmetric numerical relativity. 
That implementation has also been applied to a study of collapse of rotating supermassive stars to supermassive black 
holes [22]. 

We have recently remade a hydrodynamic implementation using a high-resolution shock-capturing scheme based on 

a Godimov-typc scheme [23 28]. Although the previous one [21] works well for problems in which shocks arc weak, 
such as evolution of single rotating stars and collapse of neutron stars and supermassive stars to a black hole, it is 
expected that such implementation cannot produce an accurate numerical result for problems in which shocks are 
strong. During rotating core collapses to a neutron star or a black hole, strong shocks are likely to be accompanied. 
Therefore, implementing a high-resolution shock-capturing scheme, such as that adopted in [28], is a promising 
strategy. To check that the new implementation works well, we have performed a wide variety of test simulations. 
In this paper, we present the numerical results, paying particular attention to long-term numerical simulations of 
neutron stars as done in, e.g., [18,28], and to stellar collapse in which strong shocks are accompanied. Finally, we 
present the first numerical results of stellar core collapse of a neutron star for which the simulation is started from 
a realistic initial condition. In addition to presenting the successful numerical results, we address the advantage of 
axisymmetric simulations in testing a new general relativistic hydrodynamic implementation, since we can study in 
detail the convergence of numerical results in the test simulations changing the grid number for a wide range to a 
well-resolved level (e.g., N ~ several hundreds), which is still difficult in three-dimensional simulations because of 
restricted computational resources. 

The paper is organized as follows. In Sec. II, we describe the formulation that we adopt. In Sec. Ill, we define global 
quantities of the system and describe the calibration method for the numerical results. In Sec. IV, we present the 
numerical results. Section V is devoted to a summary and discussion. Throughout this paper, we use the geometrical 
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units in which G = 1 = c, where G and c denote the gravitational constant and speed of hght. We use Cartesian 
coordinates, x'^ = {x,y,z), as the spatial coordinates, with r = ^y + y"^ + z"^, w = \/ x"^ + y"^, and if = ian~^[y/x). 
t denotes the coordinate time. Greek indices /j,,!/, - ■■ denote x, y, z, and t, small Latin indices i, j, - ■ ■ denote x, y, and 
z, and capital Latin indices A,B,--- denote x and z. 



II. FORMULATION 



A. Formulation for the Einstein equation 

The Einstein equation is solved in the (3+1) formulation, in which the line element is written in the form 

= (_a2 + l3kP^)df + 2pkdx''dt + -yijdx'dx^ , (2) 

where a, (3'^, and 7,j are the lapse function, shift vector, and three-metric. The three-metric is defined by 

= 9niy + n^Hu, (3) 

where is a timelike, unit-normal vector that is orthogonal to a spacelike hypersurface, and its components are 
written as (l/a, —(3'^ /a). The extrinsic curvature is defined as 

Kij = -\^nlij = ('^t^y - ^il^j - ^jf^ij ' (4) 

where t£„ is the Lie derivative with respect to n^, and Di is the covariant derivative with respect to "fjk- 

The Einstein field equations are solved using the same formulation, gauge conditions, and outer boundary condi- 
tions as in previous papers [17-20,29,30,21]: We adopt the so-called Nakamura-Shibata formulation [2,31] with some 
modification from the original version (see [20] , to which the reader may refer for basic equations and gauge conditions 
in the latest version). In this formalism, we evolve the following geometric variables using a free evolution code: 

(5) 

(6) 
(7) 

(8) 
(9) 

The Hamiltonian and momentum constraint equations are solved at f = 0, and used to check the accuracy of numerical 
solutions during computation. 

The slicing and spatial gauge conditions for determining a and (3'^ are basically the same as those adopted in our 

previous scries of papers [17-20,29,30,21], i.e., we impose an "approximate" maximal slicing condition (K w 0) and 
an "approximate" minimum distortion (AMD) gauge condition [Di{dtj^^) ~ 0, where Di is the covariant derivative 
with respect to ^ij]. However, in contrast with previous papers [17,29,19,20], we do not modify the spatial gauge 
condition even in the formation of black holes, since in the axisymmetric simulation, a sufficient number of grid points 
can be taken to resolve black hole formation and subsequent evolution even using the AMD gauge condition without 
modification. 

We impose the axially symmetric condition to the geometric variables using the so-called cartoon method proposed 
by Alcubierre et al. [16]. First, we define the computational domain as < a;,z < L and —Ay < y < Ay, where L 
denotes the location of the outer boundaries, and reflection symmetry with respect to the z = plane is assumed. 
With this computational domain, we need only three points in the y direction, and ±Ay. We determine here that 
the Einstein equation is solved only in the y = plane. Then, the boundary conditions at y = ±Ay that are necessary 
in evaluating y derivatives are supplied from the assumption of axial symmetry as 

Qab= AfAi>Q<§l, 

Q..= Qi?, Qz = gf, Q = Q^°\ (10) 



4> = 


^ln[det(7y)]. 


lij = 




K= 
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where 



^ ' sin</5(x) cosi^(a;j ' ^ ' 



and </?(a;) = tan ^[±Ay / ^/x^+{Ay)^]. Qij, Qi, and Q denote (7^, ^y), {Fi,f3^), and ((/f), ii', a) , respectively, and 
Q.-^' , Q,-"' and Q^"' arc the vahics of Qij.Qi, and Q at (-^/.x^ + (Ay)^, 0, z), which are interpolated using Lagrange's 
formula [32] with three nearby grid points along the x direction (i.e., x ± Ax and x). At x — we use only two 
points, X — Ax and x, for the extrapolation. 

To impose the gauge conditions, as well as to solve the constraint equations in preparing the initial conditions, we 
solve scalar and vector elliptic-type equations of the form [18] 

AflatQ = S, (12) 
AflatQi = Si, (13) 

where Aflat denotes the Laplacian in the flat three-dimensional space, and S and Si denote the source terms. Using 
the interpolation mentioned above, dyyQ and dyyQi are evaluated in the finite differencing as 

Oyy^-^ (A2/)2 ' OyyUl,-Z ^^^^^ 

_ ^ gi°^|cosy(3;)| - Q:,{x,Q,z) _ ^ Q^°^| cosy(a;)| - Qy{x,Q,z) 

Oyv^x - ^ (Ay)2 ' yy^y ~ {AyY 

On the other hand, the finite differencing in the x and z directions, dxxQi and d^zQi, is written in the standard form 
as 

Qi{x + Ax, 0, z) - 2Q,{x, 0, z) +Qi{x- Ax, 0, z) 
(Ax)2 

Qi{x, 0,z + Az) - 2Qi{x, 0, z) + Qi{x, Q,z- Az) 
{AzY 

Thus, in the finite differencing form for each component of Eq. (13), only one component of Qi is included, implying 
that each component of the vector elliptic-type equation is solved independently, as in the case of the scalar elliptic 
equation. 

Finally, we note a necessary modification for numerically handling Einstein's evolution equations in the axisymmetric 
case. In our formalism, the evolution equations are written in the form [31,17,20] 

dtQ + l3''dkQ = right - hand side, (14) 

where Q denotes one of the geometric variables jij, Aij, cj), K, and -Ffc. In the three-dimensional case, we apply an 
upwind scheme to numerically handle the transport term (i^dkQ for all the components [17]. In the axisymmetric 
case, the same method is used for P'^dxQ and P^dzQ, but is not for pydyQ, since it is not appropriate. (Remember 
that in the hydrodynamic equations in the axisymmetric case, there is no transport term for the rotational direction.) 
For pydyQ, we use the following schemes: For Q = cj), K, 7^^, and A^z, we set (3ydyQ = because of symmetry. For 
other variables, we simply use the cell-centered (second-order) finite difference. 

B. Formulation for the hydrodynamic equations in general relativity 
The hydrodynamic equations in general relativity are written as 

v^ipun = 0, (15) 

V^Tt = 0, (16) 

where is the covariant derivative with respect to the spacetime metric gfj,^,, p is the baryon rest-mass density, u** 
is the four-velocity, and 
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T'"' = phuV + Pg"". (17) 

Here, P is the pressure, h = 1 + e + P/pis the enthalpy, and e is the specific internal energy. Equations (15) and (16) 
are the continuity equation and the equations of motion, respectively. 

We adopt the so-called high-resolution shock-capturing scheme in numerically handling the transport terms of 
hydrodynamic equations. To use such a scheme, the hydrodynamic equations should be of a conservative form as 

dt{p*Vv) + di{p*^v^) = 0, (18) 
dt{p*^Uj) + di{p*^v% + Pae^'^^S'j) 

= Pdj{ae^'>'^) - p*^[whdja - Uidj(3' + ^yff^kUidjl''^] (19) 

dt{p*e^) + di[p*^ev' + Pe^'>'^{v' + /?')] 

= ae^'t'^PK + ^^Uiii^K'^ - p.^u^^Wja, (20) 

where 

p* = frwe^'*', (21) 

v' = ^ = -/3^ + aY'^, (22) 
w nw 

Ui = hui, (23) 

p60 p 

e = — Ti,^nV = hw , (24) 

p* pw 

w = au*, (25) 

and 77 is a determinant in curvilinear coordinates; in the cylindrical coordinates, rj = w. We note that subscripts i, j, - ■ ■ 
here denote the components in curvilinear spatial coordinates. Equations (18), (19), and (20) are the continuity, Euler, 
and energy equations. The Euler and energy equations are derived from 7"^^ V^T'^ = and n^'V^T'^ = 0, respectively. 

We solve the hydrodynamic equations in the assumption of axial symmetry. Thus, we first write equations in the 
cylindrical coordinates (tu, z). However, the Einstein equations are solved in the y = plane with the Cartesian 
coordinates. Hence, we rewrite the hydrodynamic equations in the Cartesian coordinates using relations such as 
tu = X and for y = 0. Then, the explicit forms of the equations can be written as 

dtP* + d,{p*v-) + d,{p,v'-) = -^, (26) 

X 



.P^+P^S^^ + PdAae^'^) 

(27) 



whOAOi — UjOAP H — -. — -Oat ~ 



2wh w 

dt{p*Uy) + dx{p*UyV'^) + dz{p*UyV^) = — — , (28) 



X 

dt{p*e) + 9x[p*ei;" + Pe6'^(i;" + /?")] + d,[p*ev^ + Pe^'^{v^ + /J^)] 

= _ P*^-^ + P-''(-'+n ^ „,6,p^ ^ ^u.u,K^^ - p*ua-Dja, (29) 

X u*h 

where a subscript A denotes x or z, and ■ ■ ■ here denote x,y, and z. For numerically handling the transport 
terms as dx{- ■ •) and dz{- ■ •), we apply an approximate Riemann solver with third-order (piecewise parabolic) spatial 
interpolation. Other terms are regarded as the source terms. No artificial viscosity is added, in contrast with our 
previous axisymmetric implementation [21]. The time integration is done with the second-order Runge-Kutta method 
as explained in [18]. Detailed numerical methods with respect to the treatment of the transport terms are also 
described in Appendix A. 

We note that from Eqs. (26) and (28) , conservation of baryon rest-mass and angular momentum is derived. However, 

we write these cqiiations as nonconservative forms and, hence, these conserved quantities are not precisely conserved 
in numerical computation. To suppress the growth of violation of the conservations in an acceptable level (e.g., within 
1%), we should be careful in the grid resolution (see Sec. IV). 
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In every time step of computation, w at each grid point is obtained by solving the following equation which is 
derived from the normalization relation of the four-velocity as 

w'^ = l + ^'^UiUj = 1 + ^'^UiUj (- + . (30) 

Here, P = P{p,e) = P[p>,/ {we^'^),e] (see Sec. II C) and p = p^/{we^'^). Thus for a given ^ij, </>, Ui, e, and p*, Eq. 
(30) constitutes an algebraic equation for w, which can be solved by standard numerical techniques [32] . After w is 
obtained, p, P, e, h, and can be updated. We note that this procedure is essentially the same as that used in our 
previous papers (see [18] for details). 



C. Equations of state 

We adopt two equations of state. One is the so-called F-law equation of state of the form 

P = (r - l)pe, (31) 

where T is an adiabatic constant. In using Eq. (31), we always give an initial condition using the polytropic equation 
of state P = Kpp^ of the identical F, where Kp is a polytropic constant. In this paper, we set the adiabatic constant 
as r = 2 [i.e., the polytropic index n is given by n = l/(r — 1) = 1] as a qualitative approximation of moderately 
stiff equations of state for neutron stars. We note that if we prepare a polytropic star as an initial condition and the 
system evolves in an adiabatic manner (with no shock, cooling, and heating), the equation of state is preserved in the 
polytropic form even using Eq. (31); i.e., the value P/p^ = Kp{x'^) for any fluid element remains a constant (= Kp). 

The second one is a parametric equation of state that has been used by Yamada and Sato [11], and by Miiller and 
his collaborators [12,14] for the simulation of a rotating stellar core collapse. In this equation of state, we assume that 
the pressure consists of the sum of polytropic and thermal parts as 

P = Pp + Pth. (32) 

The polytropic part is in general given as = Kp{p)p^^P\ where Kp and F are not constants but functions of 

density p. In this paper, we follow [14] for the choice of Kp{p) and T{p): For density smaller than the nuclear density 
Pnuc = 2 X lO^'^ g/cm^, F = ri(=const) is set to be < 4/3, and for p > p^uc, T = F2(= const) > 2. Thus, 



Kip^\ p<p, 
K2P^', P>P, 



,r. (33) 



where Ki and K2 are constants. Since Pp should be continuous, we demand that the relation, K2 = Kip^^^^^ , should 
be satisfied. Following [12,14], we set Ki = 5 x 10^^ cgs, because we can well approximate the polytropic part of the 
equation of state for p < p„uc in which the degenerate pressure of electrons is dominant. Taking into account that the 
specific internal energy should also be continuous at p = Pnuc, the polytropic specific internal energy ep is written as 



ep 



Ki r 

7P \ P< Pn 



^'-^ r,-i (34) 



I Fa-l^ + (ri-l)(F2-l) ' ^-^""^ 



With these settings, we mimic a realistic equation of state for high-density, cold nuclear matter. 

The thermal part of pressure plays a role in the case that shocks are generated. Here, we write it as 

Pth = (Ft - l)p£th, (35) 

where eth = e — ep- Following [12,14], we set Fth = 1-5 in this paper. 

We performed simulations of rotating stellar collapses using this parametric equation of state. In choosing this, we 
always give equilibrium stars as initial conditions using the polytropic equation of state 

P = Ko//^ (36) 

where Kq is a constant. Following [12,14], we set i^o = 5 x 10^^ cm^/s^/gr^/^, with which a soft equation of 
state governed by the electron degenerate pressure is well approximated [33]. Here, Kq and Ki are related by 
Ki = Kopy^~^^, where po = l g/cw?. 



6 



D. Adding atmosphere 

In using high-resolution shock-capturing schemes, we have to add an atmosphere of small density outside stars, 
since p and P have to be nonzero. At t = 0, wc put an atmosphere of uniform density and specific internal energy in 
the computational domain of p = 0, according to the following methods. For the F-law equation of state with F = 2, 
the uniform density of the atmosphere is set as pa = 10~^/9inax, where Pmax denotes the maximum density of a star. 
The specific internal energy is given using the polytropic constant as Kp/A. For parametric equations of state (32), 
the uniform density of the atmosphere is set as pa ~ 1 g/cm^. In this case, the specific internal energy is given using 
the polytropic constant as Kq. 

For the F-law equation of state with F = 2, the density decreases steeply around the surface of a neutron star. In 
such a case, numerical instability could often turn on around the stellar surface, if the density of the atmosphere is 
too low. This is the reason that we attach the atmosphere of relatively high density. On the other hand, a small value 
of Pa is acceptable in parametric equations of state. 



III. GLOBAL QUANTITIES AND METHOD FOR CALIBRATION 



Wc monitor the conservation of the total baryon rest-mass M*, ADM mass M, and angular momentum J, which 
are computed in the y = plane as 



M* = An 



dzp^, 



M 



J = A-K 



/ xdx 

Jo Jo 

[ p't> ^ P^4' f _ _ 2 -1 

-2 / xdx / dz -2TTEe"^ -h —R - —{AijA'^ - -K^ \ 
Jo Jo L 8 8 1. 3 J 

/ x'^dx 
Jo Jo 



dzp<tU 



(37) 
(38) 

(39) 



where E = ph/uP'—P and R is the Ricci scalar with respect to 7^ . should be conserved in any system. Because of the 
axial symmetry, J should also be conserved. On the other hand, M is not conserved in general because of gravitational 
radiation. However, the total radiated energy of gravitational waves is likely to be quite small in the axisymmetric 
spacetimc, so that wc can consider M as an approximately conserved quantity. In our axisymmetric hydrodynamic 
implementation, and J are not guaranteed to be conserved precisely. Thus, monitoring the conservation of them 
is a good check of numerical accuracy. 

In addition to the mass and angular momentum, we also check the conservation of the specific angular momentum 
spectrum [37], 



M*(jo) = An xdxdzp*, 

Jj<jo 



(40) 



where j is the specific angular momentum computed as xuy{= hu^p) and jo denotes a particular value for j. 

The numerical accuracy is also checked monitoring the violation of the Hamiltonian constraint, which is written as 



(41) 



where 'ip — e*^, and A denotes the Laplacian with respect to ^ij. In this paper, we define the averaged violation 
according to 



ERROR : 



1 



p*\V\d^x, 



(42) 



where 



V 



o o 



12' 



(43) 



+ 2nEi,^ + ^AijA'^ + 
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Namely, wc use as the weight factor for the average. The reason that we introduce this average is as follows. In 
using high-resolution shock-capturing schemes, we add an atmosphere of small density outside neutron stars and/or 
collapsing stars. In the atmosphere, a small error in the metric results in a large violation of the Hamiltonian 
constraint because is a very small value. Furthermore, the volume fraction occupied by the atmosphere in the 
whole computational domain is larger than that for main bodies. Thus, if we simply compute the volume integral of 
\V\, it is close to unity irrespective of the grid resolution. However, the numerical accuracy in the atmosphere is not 
very important for evolution of the main bodies and for global evolution of the system in which we are interested. 
Therefore, to monitor whether the main bodies (neutron stars and collapsing stars) are accurately computed or not, 
this type of weight factor is necessary. 

IV. NUMERICAL RESULTS 

In the numerical simulations reported in Sees. IV A-C below, we adopted a fixed uniform grid, in which the grid 

spacing Ax = Ay — lS.z is constant, with grid size (A^ + 1,7V + 1) for (x, z) to cover a computational domain as 
< X, z < L, where L = NAx. In the simulation reported in Sec. IV D, we varied the grid spacing during the 
computation, but still used the uniform grid in which Ax = Ay = Az. To check the convergence of the numerical 
results for Ax — > 0, numerical computations were carried out with three levels of the grid resolution while fixing 
L. All the computations were done on the FACOM VPP5000 machine in the data processing center of the National 
Astronomical Observatory of Japan. The memory and CPU time in one run with a grid size N = 600 and with four 
processors axe about 1 GByte and 12 CPU hours for ~ 30000 time steps. 

We note that for a simulation with N = 600 in three spatial dimensions, we would need ~ 300 Gbytes memory and 
it takes ^ 1000 CPU hours for 30000 time steps using 32 processors [20]. Such simulation is pragmatically impossible, 
because the computational costs are too high (note that under normal circumstances, we can use at most ~ 1000 
CPU hours per year). However, taking N = 600 in an axisymmetric simulation is an easy task with the current 
computational resources. 

A. Spherical neutron stars 

In this subsection, we focus on the long-term evolution of spherical neutron stars. The initial conditions are given 
using the polytropic equation of state with F = 2, and during the time evolution, the F-law equation of state (31) 
is used. Although we did the same test simulations for the previous hydrodynamic implementation and obtained 
successful outputs from it [18], we repeated the tests again in the present new implementation to demonstrate that it 
also works well. A difi^erence in the present tests from the previous ones is that we have performed much longer-term 
simulations than those in the previous tests, since it is computationally inexpensive and pragmatically possible to do 
in the axisymmetric case. 

In the polytropic equations of state, the polynomial relation c^~"ifp^^G~^/^ has dimension of mass. With this 
property, all the quantities can be scaled to be nondimensional, if we multiply an appropriate combination of c, G, 
and Kp. Thus, we will only show the nondimensional quantities in using this equation of state. In other words, we 
adopt the units of c = G = Kp = 1. In these units, the maximum ADM mass and baryon rest-mass of spherical 
neutron stars are w 0.164 and 0.180, respectively, with the central density pc ~ 0.318. Recovering Kp, the mass and 
density in dimensional units can be written as 

/ Kp V'^ f \ 

...„ = 1.35xlO-g/cn.3(^_i|_)-'(iL) 

In Table I, we list several quantities for five models of spherical neutron stars that we pick up in this subsection. 
Below, we refer to these models as models (SI) - (S5). Models (SI) (S4) are stable against gravitational collapse, 
while (S5) is marginally stable. 

In Figs. 1 (a) and (b), we display the time evolution of the central density, central value of the lapse function 
(hereafter ac)-, central value of K (Kc), averaged violation of the Hamiltonian constraint, and violation of ADM mass 

— 1/2 

and baryon rest- mass conservation for model (S2). Throughout this subsection, the time is shown in units of Pcinit' 
where Pc.init denotes the central density at t = 0. To induce a small oscillation, we initially reduce the pressure 
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(44) 
(45) 





Pc 




M 


M/R 


P o"'^ 


SI 


0.0637 


0.105 


0.100 


0.0932 




S2 


0.127 


0.150 


0.140 


0.146 


5.0 


S3 


0.191 


0.170 


0.156 


0.178 


6.9 


S4 


0.255 


0.178 


0.162 


0.200 


11 


S'>' 


0.318 


O.ISO 


O.Kil 


0.21 1 





TABLE I. The central density pc, baryon rest-mass Af, , ADM mass M, and compactness M/R of spherical neutron stars 
with r = 2 that we pick up in this paper. Here, R denotes the circumference radius. All the quantities are shown in units 
of c = G = Kv = 1. The star denoted with f mark is of the maximum allowed mass and hence marginally stable against 
gravitational collapse. In the last column, numerical results of the radial oscillation period for the / mode. Pose, are presented 
in units of Pc^^^ . 
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(a) (b) 
FIG. 1. (a) Time evolution of central density, central value of lapse function, and extrinsic curvature at origin, and (b) time 
evolution of averaged violation of the Haniiltonian constraint, violation of the ADM mass conservation, and violation of the 
baryon rest-mass conservation for a stable spherical neutron star (S2). In both figures, the solid, dotted, and dashed curves 
denote the results with N = 180, 120, and 90. With these grid numbers, the radius is initially covered by about 65, 44, and 33 
grid numbers. Pc.init denotes the central density at t = 0, and the time is shown in units of pjf^t- 



by 0.2%. We note that whenever we superimpose a perturbation to an equihbrium configuration, we reinforce the 

Hamiltonian and momentum constraints at t = 0. Numerical results are shown for A'^ = 90, 120, and 180 with a 
fixed value of L, to demonstrate that the convergence is achieved. The simiilations continued for ~ 30 dynamical 
time scales (see also Fig. 2) until the crash of the run, irrespective of the grid resolution, although the accuracy 
deteriorates gradually with time. Here, we refer to the period of the fundamental radial (and quasiradial) oscillation 
as the dynamical time scale. 

L is chosen as ~ SRg, where Rs denotes the coordinate radius of the neutron star. For simulating spherical systems, 
we imposed the outer boundary conditions as 

% = 6ij, Aij = 0, {r(l)),r = 0, K = 0, and Fi = 0. (46) 

(For nonsphcrical problems, wo impose an outgoing boundary condition for jij and Aij.) These boundary conditions 
are adequate but not physically perfect. As a result, for the choice of a too small value of X as ~ the numerical 
solution is affected by the spurious effects of the outer boundaries, resulting in an earlier crash of the run. However, 
for L > 2Rs, the results are not significantly modified by the spurious effects. 

As we mentioned in Sec. Ill, the baryon rest-mass is not numerically conserved strictly. However, the violation 
does not seriously affect the numerical results. Indeed, the averaged values of the central density and lapse remain 
constants, as they should. The numerical results indicate that if we want to suppress the violation of the mass 
conservation within 1% (2%) after 10 dynamical time scales, the radius of the neutron star should be covered by more 
than 40 (30) grid points. 
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(a) (b) 
FIG. 2. (a) Time evolution of central density, lapse function at origin, and averaged violation of the Hamiltonian constraint 
for stable stars (S2) (solid curves), (S3) (dotted curves), and (S4) (dashed curves). The simulations were performed with 
N = 180. With these grid numbers, the radius is covered by about 65, 58, and 51 grid numbers for (S2), (S3), and (S4), 
respectively, (b) Fourier spectra of the central density for (S2) (solid curve), (S3) (dotted curve), and (S4) (dashed curve). 



The error of mass conservation converges to zero with improving the grid resolution at approximately second- 
order. The averaged violation of the Hamiltonian constraint also indicates approximate second-order convergence. 
Therefore, we can conclude that the numerical solution converges to the exact solution in the limit Ax 0. We note, 
however, that the convergence is only approximately at second-order, because near stellar surfaces, the gradients of 
hydrodynamic variables are so steep that transport terms arc often computed with first-order accuracy in space. The 
convergence may also become first-order if shocks are generated during numerical computations (see Sec. IV C), since 
near the shocks, the hydrodynamic computations are done with first-order accuracy. A similar tendency is reported 
by Miller et al. in the simulations of a head-on collision of two neutron stars [34] . 

In the last figure of Fig. 1(a), it is shown that Kc is not zero exactly but relaxes to a finite value. This indicates 
that even in solving the equation for the maximal slicing condition, K deviates from zero as long as finite differencing 
methods are used. Indeed, Kc converges to zero at second-order with improving the grid resolution. Thus the maximal 
slicing condition K = Q cannot be precisely imposed in numerical computation even in a well-resolved simulation, if 
we adopt finite difi^erencing schemes. This tells us that we should follow the evolution of K and should not a priori 
set K = Q la numerical computation, even in choosing the maximal slicing condition. 

Long-term simulations for more compact stable stars (S3) and (S4) were also carried out. In Fig. 2, we display the 
results for models (S2), (S3), and (S4) together. For (S3), the simulation continued for ^ 20 dynamical time scales 
until the crash of the run. However, for (S4), the star starts collapsing to a black hole after about five oscillation 
periods. The reason for this consequence is clear. The ADM mass of model (S4) is « 99% of the maximum allowed 
value. Thus, with a slight increase of the mass as a result of the accumulation of numerical error, the mass exceeded 
the maximum allowed value for stable stars, resulting in eventual gravitational collapse. It is interesting to note that 
in this case, the computation was able to be continued until the formation of a black hole of the apparent horizon mass 
[see Eq. (49) for definition] ^ M, where M is the initial ADM mass. Thus, the increase of the central density and 
the decrease of etc [see Fig. 2 (a)] do not imply that the computation crashed. To avoid the collapse to a black hole 
and to make the oscillating time longer, we need to take more grid numbers to improve the grid resolution. Actually, 
we have checked that we can increase the oscillation cycles in numerical computation with improvement of the grid 
resolution. 

The duration of the simulation to crash for model (S3) was shorter than that of (S2). This fact indicates that with 
increasing the compactness of neutron stars, the computations crash earlier. We note here that the duration of the 
simulation for these models depends only weakly on the grid resolution. Thus, it does not depend on the numerical 
accuracy, but seems to depend on certain factors associated with formulation or gauge conditions or outer boundaries. 
The same tendency is found in the simulations of rotating neutron stars. We will argue this point in Sec. IV B again. 

Although the simulations can be continued only for a finite time scale, the duration of > 10 dynamical time scales 
seems to be sufficiently long. Indeed, we can accurately extract the oscillation frequencies of the fundamental radial 
mode from these simulations. In Fig. 2(b), we display the Fourier spectra of Pc(t), which is defined as 
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(a) (b) 
FIG. 3. (a) Time evolution of central density, lapse function at origin, and violation of baryon mass conservation, and (b) 
time evolution of mass of the apparent horizon Mah in units of the ADM mass of the system for a marginally stable spherical 
neutron star (S5). To destabilize, we initially reduce the pressure by 0.5%. In both figures, the solid, dotted, and dashed curves 
denote the results with N = 180, 120, and 90. With these grid numbers, the radius is initially covered by about 71, 47, and 35 
grid points. 



F(/)= /'bc(i)-Pc,av]e'"'^*rft , (47) 
Jo 

3^/2 1/2 1/2 

where tf is chosen as ~ 70pc , 55pc , and 40pc for (S2), (S3), and (S4), respectively. pc,av is computed from 

if Jo 

2/2 1/2 

The Fourier spectra indicate that the oscillation period of the / mode of the radial oscillation is ~ h.Qpc , 6.9pc , 

1 /2 

and llp^ ' for (S2), (S3), and (S4). These values agree well with those derived from Chandrasekhar's semianalytic 
formula [35,18]. Furthermore, the result for (S2) is in good agreement with that for a spherical star of = 0.128 
reported in [28]. Thus, we conclude that the computation can be continued for a sufficiently long time to accurately 
obtain the oscillation frequencies of even extremely relativistic neutron stars. The simulation would be able to be 
carried out to study nonspherical oscillations of neutron stars, as was done in [18,36]. 

In Fig. 3, we display the time evolution of several quantities for collapse of a marginally stable spherical neutron 
star (85). To induce the collapse, we initially reduced the pressure by 0.5%. With this setting, the neutron star 

— 1/2 

collapses to a black hole in ~ lOp^ jj;;^. We have checked that even with a 0.2% decrease of the pressure, the star 

— 1/2 

collapses to a black hole in ~ 20p^ which is longer than that for the case of 0.5%. Numerical results are presented 
for N — 90, 120, and 180 and demonstrate that the convergence is achieved. Note that with lower grid resolution, it 
takes a slightly longer time to form a black hole. This is because the dissipation that prevents the increase of density 
is larger with lower grid resolution. 

In the final phase of the collapse, the grid resolution around the black hole forming region became so bad that 
the computation crashed. It is interesting to note that the magnitude of the averaged violation of the Hamiltonian 

— 1/2 

constraint, ERROR, relaxes to ~ 0.2 irrespective of the grid resolution at t ~ lOp^i^^. This verifies that the 
computation crashes when ERROR is ~ 0.2. We can also observe that as the accuracy deteriorates, (i) the central 
density stops increasing and instead starts decreasing, and (ii) exponential decrease of ac that is a feature in the 
maximal slicing condition is modified. The reason for (i) is as follows. In our implementation, is a fundamental 
quantity to evolve, and p is computed from p^/{e*"'''w). In numerical computations, and (f) monotonically increase, 
but since 4> around the origin is too large [of 0(1)] in the late phase of the collapse, a small error in (j) leads to a large 
error in p. As a result, p decreases in the late phase. The reason for (ii) is simply that the computation crashed. 
Indeed, the time at which the behavior of ac starts changing agrees with that at which the magnitude of the averaged 
violation of the Hamiltonian constraint saturates to ^ 0.2. 

Although the accuracy deteriorates in the final phase of the collapse, the simulation can be carried out at least until 
the formation of an almost static black hole. To confirm the black hole formation, apparent horizons were located 
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during the simulations. In Fig. 3(b), we display the time evolution of the mass of the apparent horizon in units of 
the ADM mass of the system. Here, the mass of the apparent horizon is defined as 



M.h = Vt6^' (49) 

where S is the area of the apparent horizon. The figure indicates that Mah relaxes approximately to M in the final 
phase of the collapse. (For N = 180, \Mau/M — 1| is less than 1%.) This implies that the simulation was carried out 
up to the time when a spacetime settles down to a static black hole spacetime. 




FIG. 4. (a) Time evolution of central density, lapse function at origin, and averaged violation of Hamiltonian constraint for 
an oscillating spherical neutron star. The solid, dotted, and dashed curves denote the results with A' = 180, 120, and 90. With 
these grid numbers, the radius is initially covered by about 78, 52, and 39 grid points, (b) Profiles of e/{Kpp^~^) and p at 
t w (dotted curves) and after about one oscillation period (solid curves). The unit of the horizontal axis is the initial radius 
of the neutron star. 



In Fig. 4(a), we show the time evolution of several quantities for an oscillating spherical neutron star of a high 
amplitude. In this simulation, we picked up a low-mass spherical star (SI), and to induce an oscillation of a high 

amplitude, wc initially reduced the pressure by 40%. The numerical results are shown for N = 90, 120, and 180, and 
demonstrate that convergence is achieved. In each case, the neutron star is initially covered by 39, 52, and 78 grid 
points, respectively. 

Because of a significant decrease of the pressure, the radius of the neutron star decreases by a factor of 2 soon after 
the simulation starts. However, the magnitude of the pressure decrease is not large enough for the star to collapse to 
a black hole. Instead, the star bounces when the central density becomes about six times of initial value, and repeats 
oscillations subsequently. As the density approaches the maximum, shocks are formed around the stellar surface, and 
as a result, outer envelopes explode. To illustrate that shock heating indeed occurs, we display £/{Kpp^~^) at t « 
and after about one oscillation period in Fig. 4(b). As mentioned in Sec. II C, in the absence of shocks, this quantity 
does not change from the initial value, but in the presence of shock heating it increases. Figure 4(b) clearly shows 
that in the outer envelope, the shock heating is significant. On the other hand, a negligible effect of shocks can be 
seen around the central region. Since the shocks arc generated only in the atmosphere, the averaged violation of the 
Hamiltonian constraint still converges approximately at second-order with improving the grid resolution. 

Figure 4(a) indicates that the amplitude of the oscillation gradually decreases, and after several oscillation periods, 
it settles down approximately to a constant. This illustrates that the kinetic energy of the oscillation is dissipated 
by the shocks gradually. Similar results are reported in [28] for a simulation of a migrating neutron star. The lower 
figure of Fig. 4(b) shows the density profiles at t = and t ~ one oscillation period. This indicates that the density 
profile is modified to a more centrally condensed state as a result of the shock dissipation. 

We emphasize that this simulation is highly dynamical and general relativistic. Even in such a case, the simulation 
was able to be continued for more than 20 oscillation periods. This illustrates the robustness of our implementation 
for dynamical problems in general relativity. 
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Pc 


M* 


M 


Mj ti 




JIM 


\TjW\ 


^osc Pc 


RO 


0.103 


0.169 


0.158 


0.111 


8.52 


0.667 


0.0932 


5.7 


Rl 


0.136 


0.186 


0.172 


0.129 


8.53 


0.630 


0.0909 


6.5 


R2 


0.183 


0.200 


0.183 


0.148 


8.56 


0.599 


0.0876 


8.2 


R3 


0.215 


0.204 


0.186 


0.158 


8.60 


0.585 


0.0856 


10 


R4 


0.253 


0.206 


0.188 


0.167 


8.65 


0.573 


0.0834 


16 


R5+ 


0.296 


0.206 


0.188 


0.175 


8.72 


0.561 


0.0809 





TABLE II. The central density, baryon rest-mass, ADM mass, M/R, rotational period (Prot) in units of p'l^i^, J/M^, and 
|T/W| of rotating neutron stars at mass shedding limits with F = 2 that we pick up in this paper. Here T and W are the 
rotational kinetic energy and gravitational potential energy, and R denotes the circumference radius at the equator. All the 
quantities are shown in units of c = G = Kp = 1. The star denoted with t is unstable. In the last column, numerical results 
of the quasiradial oscillation period for the / mode. Pose, are presented in units of p^ 



B. Rapidly rotating neutron stars 



We focus here on the long-term evolution of rigidly and rapidly rotating neutron stars at mass shedding limits for 
which the angular velocity at the equator is equal to the Kepler angular velocity. Following Sec. IV A, we adopt 
the polytropic equation of state with F = 2 for setting initial conditions and evolve neutron stars using the F-law 
equation of state (31). The axial ratio of the polar radius to the equatorial one is « 0.58 for rotating neutron stars 
at mass shedding limits with F = 2. As in Sec. IV A, we only present the scaled dimensionless quantities with the 
units of c = G = Kp = 1 throughout this subsection. In these units, the maximum ADM mass and baryon rest-mass 
of rigidly rotating neutron stars are about 0.188 and 0.207, respectively, with the central density pc ~ 0.27 [38]. The 
central density of the marginally stable star against gravitational collapse has slightly larger density (w 0.295) than 
this value [38,29]. In Table II, we list several quantities of six rotating neutron stars that we pick up here. In the 
following, we refer to these neutron stars as models (RO) - (R5). Models (RO) - (R4) are stable against gravitational 
collapse, and (R5) is unstable and very close to the marginally stable point. 
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(a) (b) 
FIG. 5. (a) Time evolution of central density, lapse function at origin, and extrinsic curvature at origin, and (b) time evolution 
of averaged violation of the Hamiltonian constraint, violation of rest-mass conservation, and violation of angular momentum 
conservation for a stable and rapidly rotating neutron star at the mass shedding limit (Rl). In both figures, the solid, dotted, 
dashed, and dotted-dashed curves denote the results with N = 240, 180, 120, and 90. With these grid numbers, the polar 
(equatorial) radius is covered by about 46 (80), 35 (60), 23 (40), and 18 (30) grid points, respectively. 



In Fig. 5, we display the time evolution of several quantities for model (Rl). The ADM mass of this model is 
« 91% of the maximum mass of rigidly rotating neutron stars with F = 2, so that it is a sufficiently relativistic 
model. To induce a small oscillation of the fundamental quasiradial mode, we initially reduce the pressure by 0.5%. 
Numerical results are shown for N = 240, 180, 120, and 90 to demonstrate that convergence is achieved. With these 
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(a) (b) 
FIG. 6. (a) Time evolution of central density and lapse function at origin for stable rotating stars (RO) (dotted-dashed 
curves), (Rl) (solid curves), (R2) (dotted curves), (R3) (long-dashed curves), and (R4) (dashed curves). The simulations were 
performed with N = 180. With these grid numbers, the equatorial radius is covered by about 60 grid numbers, (b) Fourier 
spectra oi pc{t) for (RO) (dotted-dashed curves), (Rl) (solid curves), (R2) (dotted curves), (R3) (long-dashed curves), and (R4) 
(dashed curves). 



grid numbers, the polar (equatorial) radius is covered by about 46 (80), 35 (60), 23 (40), and 18 (30) grid points, 
respectively. 

The simulations continued for ~ 10 dynamical time scales, and eventually crashed. The duration of the simulation 
to crash depends only weakly on the grid resolution as in the spherical cases. This implies that the crash is not 
triggered by accumulation of the numerical error. The duration also does not vary much even if we change the outer 
boundary conditions and the location of the outer boundary as L w 3i?e-4i?e, where Re denotes the coordinate radius 
at the equator. Furthermore, the duration of the simulation for a rotating neutron star is shorter than that for a 
spherical star of identical compactness. Therefore, we deduce that the crash of computations might be associated 
with our choice of the spatial gauge condition or the formulation, although we do not understand the reason fully 
at present. There may still be room to improve the spatial gauge condition and/or the formulation, if one wants to 
perform an extremely long-term simulation of the duration l^- 10 dynamical time scales. However, 10 dynamical time 
scales are long enough to produce scientific results for most problems, so that we do not address this problem any 
longer in this paper. 

As in the spherical case, convergence is achieved with improvement of the grid resolution. As argued in Sec. Ill, the 

angular momentum as well as the baryon rest-mass are not conserved strictly, although they are conserved quantities. 
However, the violation converges to zero nearly at second-order with improving the grid resolution. The results of 
this convergence test indicate that the polar axis should be covered by at least ^ 30 grid points if one wants to 
demand that the violation of the conservation of angular momentum and baryon rest-mass is less than a few % after 
~ 10 dynamical time scales. If the polar axis is covered by fewer than 20 grid points, the magnitude of the violation 
becomes larger than 10% after 10 dynamical time scales. 

The long-term simulations were also performed for models (RO), (R2), (R3), and (R4). In all these simulations, 
we initially reduced the pressure by 0.5% and take N = 180. With this grid number, the polar (equatorial) radius 
is covered by 35 (60) grid points. The time evolution of the central density and central value of the lapse function 
(ccc) are shown together in Fig. 6(a). As in the spherical case, the duration of the simulation is shorter for more 
compact stars. (The moment of the crash of a run is identified with the time at which Uc sharply drops.) This might 
be evidence that with our spatial gauge, coordinate distortion is accumulated too much for the long-term simulation, 
because it is likely to be accumulated more rapidly for more compact stars. As in the simulation for (S4), a high-mass 
rotating star (R4), for which the ADM mass is « 99.5% of the maximum, collapses to a black hole in a few dynamical 
time scales instead of crashing, because of a slight increase of the baryon rest-mass due to numerical error. It is 
necessary to take N > 180 to continue computations of such high-mass stars for more than two oscillation periods. 
However, as long as the neutron star is not very close to the marginally stable point, the simulation can be continued 
for more than five dynamical time scales with N ~ 200, and this duration is long enough to produce scientific results 
for most problems. For example, from these simulations, we can extract the frequency of fundamental quasi-radial 
oscillation modes. In Fig. 6(b), we show the Fourier spectra of the central density as in the case of Sec. IV A. The 
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figure indicates clear peaks that denote the fundamental frequency of the quasiradial oscillation modes. 

In Fig. 7, we summarize the frequencies of the radial and quasiradial oscillation modes for spherical and rapidly 
rotating neutron stars with F = 2. The filled and open circles denote the numerical results for spherical and rotating 
neutron stars, respectively. The filled circle along the vertical axis is plotted according to a Newtonian analysis for 
the spherical polytrope [39]. On the other hand, the filled and open circles along the horizontal axis are plotted by 
the fact that the frequency of the fundamental radial and quasiradial oscillations of the marginally stable stars is zero. 
We note that the frequency in dimensional units is computed from 

/ . 9485 H.(P„..;£,.,- ( j^)""^ {^) "\ m 

where Pose = 1// is the oscillation period of the hmdamental quasiradial mode. 

As mentioned in Sec. IV A, the frequencies for spherical neutron stars agree well with semianalytical results [35]. 
For rotating neutron stars, the frequency is slightly smaller than that for spherical neutron stars for identical Pc- A 
similar tendency is found in the results of [28] . It is interesting to note that for the identical value of the central lapse 
function, the frequencies approximately coincide for ac ^ 0.6. 




(a) (b) 
FIG. 7. Frequency of fundamental radial and quasiradial oscillation modes for spherical stars (filled circles) and rotating 
stars at mass shedding limits (open circles) (a) as a function of the central density and (b) as a function of the lapse function 
at origin. The frequency at Pc, init = for the spherical star is derived in a Newtonian analysis. The filled and open circles 
along the horizontal axis are plotted according to the fact that the oscillation frequency of the marginally stable star is zero. 



In Fig. 8(a), we display the central density, central value of the lapse function, and averaged violation of the 
Hamiltonian constraint for collapse of an unstable rotating neutron star (R5). To induce the collapse, we initially 
reduced the pressure by 0.5%. We also carried out the simulations with the reduced factor of 0.2%, and have found 
that the star (R5) collapses also in this case, although it takes longer to be a black hole. Numerical results are shown 
for N = 240, 180, and 120 and demonstrate that the convergence is achieved well. With these grid numbers, the polar 
(equatorial) radius is initially covered by about 69 (120), 52 (90), and 35 (60), respectively. 

During the collapse, the density (lapse function) monotonically increases (decreases) with time, and finally a black 
hole is formed (i.e., the apparent horizon is located) in the late time when ac 0.03. As in the spherical case, the 
reason that pc decreases in the late stage of the collapse is that a small error in (j) that is of 0(1) at the formation 
of the black hole leads to a large error in p that is computed from p^/{we^^). It is also found that at the time when 
the magnitude of the averaged violation of the Hamiltonian constraint becomes 0.2, the computation crashed due 
to the grid stretching around the horizon of a black hole, and that the magnitude of ERROR at the crash is almost 
independent of the grid resolution. 

In Fig. 8(b), we show the time evolution of the mass of the apparent horizon defined by Eq. (49). Here, the 
long-dashed horizontal line denotes the expected final value \/ S'eh/ (IGttM^), where ^eh is derived from the formula 
for the area of the event horizon of a Kerr black hole as [33] 

5EH = 87rM2(^l + yi-^), (51) 
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(a) (b) 
FIG. 8. (a) Time evolution of central density, lapse function at origin, and averaged violation of the Hamiltonian constraint, 
and (b) mass of the apparent horizon as a function of time for collapse of an unstable and rapidly rotating neutron star at the 
mass shedding limit (R5). In both figures, the solid, dotted, and dashed curves denote the results with N = 240, 180, and 120. 
With these grid numbers, the polar (equatorial) radius is initially covered by about 69 (120), 52 (90), and 35 (60) grid points. 



where M and J are the ADM mass and angular momentum of the collapsing neutron star. Since almost all the 
matter eventually falls into a black hole in this simulation, the area of the apparent horizon should settle down to 
Seh- The figure indicates that the area of the apparent horizon asymptotically approaches the expected value. This 
demonstrates that the spacetime in the final phase of our simulation almost relaxes to a stationary, Kerr black hole 
spacetime. 



C. Collapse of rotating stars with parametric equations of state 

The purpose of this subsection is to demonstrate that with our implementation, it is feasible to carry out stable and 

accurate simulations for the collapse of rotating stars with parametric equations of state (32) that are more realistic 
for high-density matter than the F-law equation of state used in the previous two subsections. Since the equation 
of state for high-density matter is still not precisely known, the parametric equation of state (32) is used for several 
choices of Fi and Initial conditions are set up adopting the polytropic equation of state (36) with F = 4/3. 

In the realistic core collapse of massive stars, the central density just before the collapse is of order 10^° g/cm^ 
[14,40]. Since the collapse leads to the formation of a neutron star of density of order lO'^'"' g/crn*' or a black hole, 
the characteristic length scale changes by a factor of ~ 100. This implies that we need to take of 0(10'^) for 
a well-resolved simulation in the fixed uniform grid. Although it is possible to take a large value of N as several 
thousands, performing such large-scale simulation is not computationally inexpensive even in the axisymmetric case. 
Since the main purpose of this subsection is not to present scientific results, but both to demonstrate that realistic 
equations of state can be adopted in our implementation and to grasp characteristic behaviors associated with new 
implementation with such equations of state, here we pick up more compact stars of central density « 6 x 10"'^^ g/cm^ 
as initial conditions to save the computational costs as a first step. In the next subsection, we will show a numerical 
result with a more realistic initial data set as illustration. 

Velocity profiles of equilibrium rotating stars used as initial conditions are given according to a popular relation 
[41,42], 

U'U^ = Vulina - n), (52) 

where ila denotes the angular velocity along the z axis, and zud is a constant. In the Newtonian limit, the rotational 
profile is written as 

^ = »a J^} 2 - (53) 
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TABLE III. The central density, baryon rest-mass, ADM mass, equatorial radius Re, J/M^ , and IT/M^I of initial conditions 
for the simulations of stellar collapse in Sec IV C. 







r2 




M(Me) 


M»corc(M0) 


Pc (g/cm-') 


a 


1.325 


2 


1.425 


1.363 


1.362 


2.68dl5 


b 


1.3 


2 


0.991 


0.925 


0.980 


6.01dl5 


c 


1.325 


2.5 


2.259 


2.056 


2.183 


1.69dl5 


d 


1.3 


2.5 


1.810 


1.600 


1.792 


2.87dl5 



TABLE IV. Maximum baryon rest-mass, ADM mass, core rest-mass M.core, and density at the maximum for spherical stars 
of cold, parametric equations of state (33), with four types of Fi and 



Thus, Wd indicates the steepness of differential rotation. In this paper, we pick up the rigidly rotating cases in which 

vjd OO and a differentially rotating case in which vJd/Re = 1/2, where i?,. is the equatorial coordinate radius. In 
both cases, we choose the axial ratio of polar radius to equatorial radius as 2/3. In Table III, we list several quantities 
for the models we adopt in the numerical computation. We refer to these models as (CI) and (C2). 

The simulations were carried out for (ri,r2) = (1.325,2) (a), (1.3,2) (b), (1.325,2.5) (c), and (1.3, 2.5) (d). In the 
following, we specify our choice of Fi and r2 in terms of (a), (b), (c), and (d). For example, if we pick up model 
(CI) as the initial condition and choose Fi = 1.325 and F2 = 2, we refer to this model as (Cla). In Table IV, we list 
masses and central density for spherical neutron stars of maximum mass in the parametric, cold equation of states 
(33) with four choices of Fi and which are obtained by numerically solving the TOY equation [33]. For Fi = 1.3 
and F2 = 2, the maximum mass is too small to be a realistic value, but with this model, we can study qualitative 
properties of stellar core collapses in an extremely soft equation of state. 

According to a numerical result [43], rigidly rotating stars in equilibrium with F = 4/3 are unconditionally unstable 
against gravitational collapse, if the compactness M/Rp is larger than ^ 1/700. Thus, the equilibrium state of model 
(CI) is unstable. The criterion of the instability has not been established yet for differentially rotating stars. Thus, 
the stability is not clear for (C2). However, since the compactness is much larger than 1/700, the initial equilibrium 
of (C2) is also likely to be unstable. Hence, by decreasing F from 4/3 to Fi < 4/3 at i = 0, the collapse is accelerated. 
To investigate convergence, the simulations were carried out with N = 600, 400, and 300 for all the models. In these 
simulations, the equatorial radius was initially covered by N grid points. 



t=2.80 msec t = 3.98 msec t = 5.16 msec 




X(km) X(km) X(km) 



FIG. 9. Density contour curves of p and velocity fields of v"^ for model (Cla) at selected time steps. The contour curves are 
drawn for p/pnuc = lO-'' '^-' , for j = 0, 1, 2, ■ • • , 20. 

In Figs. 9 and 10, we display snapshots of the density contour curves and velocity fields at selected time steps for 
models (Cla) and (C2a) as examples. It is found that model (CI) has a spheroidal structure initially, while model 
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FIG. 10. The same as Fig. 9 but for model (C2a). 



(C2) has a slightly toroidal shape duo to the effect of differential rotation. In Figs. 11 and 12, we also show the 
time evolution of several quantities for models (Cla) and (C2a). Here M*core denotes the total baryon rest-mass of 
high-density matter with p > pnuc- As indicated in these figures, the collapses proceed in the following manner. In the 
early phase in which the central density is smaller than /9„uc, the central region of the star contracts approximately 
in a homologous manner because the adiabatic index is close to 4/3 [44]. After the central density exceeds Pnuo a 
core is formed at the central region, and the mass gradually increases as a result of subsequent accretion. At the time 
when the central density becomes ^ 4.5pnuc for (Cla) and ^ 3.5/3„uc for (C2a), the increase of the core mass stops, 
and strong shocks rapidly propagate outward due to the restoring force of the core. After this moment, the core 
settles down toward an approximately stationary state, and finally the central density relaxes to ~ 3/9„uc for (Cla) 
and ~ 2.5/9nuc for (C2a). On the other hand, the shocks propagate outward, sweeping the infalling matter. In both 
cases (Cla) and (C2a), the baryon rest-mass of p > pnuc is ~ O.65M0 after shock propagates far from the core. Since 
the angular momentum at f = for (C2a) is larger than that for (Cla), the baryon rest-mass of the core and central 
density of the final state of the core are slightly smaller for model (C2a). These facts indicate that the products after 
collapse depend basically on the equation of state, but are modified by the magnitude and distribution of angular 
momentum that are initially retained by a precollapse star. If the initial angular momentum is much larger than those 
for (Cla) and (C2a), the collapse would be halted before the central density reaches Pnucj and as a result, a neutron 
star would not be formed [14]. 

Another noticeable feature is that for the differentially rotating initial condition, a small modulation is found in the 
time evolution of the central density and ac after a quasistatic core is formed. This is likely because in the collapse 
with a difiierentially rotating initial condition, the collapsing star deviates highly from spherical symmetry and shocks 
are formed in a nonspherical manner. As a result, the nonspherical oscillation is excited in the formed core. Similar 
results are observed in [14]. 

In Figs. 11(b) and 12(b), we display the violation of the baryon rest-mass conservation, angular momentum 
conservation, and averaged violation of the Hamiltonian constraint. First, we note that at t ~ 6.5 msec, the shocks 
reach the outer boundaries, and the matter escapes from the computational domain. This causes the rapid angular 
momentum decrease for t ^ 6.5 msec. Besides this, the errors converge to zero with improving the grid resolution, as 
we saw in Sees. IV A and IV B. However, there are several noticeable properties that have not been seen in the previous 
subsections. One is that the magnitude of the errors quickly increases when shocks are generated. The second is that 
the averaged violation of the Hamiltonian constraint converges to zero at about first-order (not at second-order) . We 
cannot explain the reason for them correctly. However, we deduce it as follows: At the shocks, the hydrodynamic 
quantities are discontinuous, and as a result, the second partial derivatives of the metric are discontinuous and the 
complete analyticity for the metric is violated. This may change the global convergence property from second-order 
to first-order. 

The magnitude of the errors for (C2a) is smaller than that for (Cla) with identical grid spacing. This is a consequence 
of the fact that the central density for (C2a) after shock formation is smaller than that for (Cla), and as a result, the 
grid resolution for (C2a) is better than that for (Cla). 

To see the effects of the equations of state on the dynamics of collapse and on accumulation of the numerical 
errors, we display numerical results for (Cla), (Clb), (Clc), and (Cld) in Fig. 13. All the simulations were done 
with N = 600, by which the equatorial radius is initially covered by N grid points. Reflecting the difference of Fi 
and the products after the collapse and time evolution of the numerical errors are different. The following is a 
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(a) (b) 
FIG. 11. (a) Time evolution of central density, lapse function at origin, and fraction of baryon rest-mass for p > pnuc, 
and (b) time evolution of violation of rest-mass conservation, angular momentum conservation, and averaged violation of the 
Hamiltonian constraint for collapse of a star (Cla). In both figures, the solid, dotted, and dashed curves denote the results 
with N = 600, 400, and 300. With these grid numbers, Aa; = 0.442, 0.663, and 0.884 km. For t > 6.5 msec, shocks reach outer 
boundaries of the computational domain, and matter starts escaping. 



summary of the differences: (a) for Fi = 1.3, the magnitude of the errors is larger than that for Fi = 1.325, (b) for 
the identical value of Fi, the central density and M^^ore in the final state of the core are smaller for the larger value 
of (c) M*core IS larger for the larger value of Fi with the identical value of (d) the magnitude of the error of 
the baryon rest-mass conservation always increases with time (never decreases in the absence of mass ejection from 
the computational domain) irrespective of Fi and F2, and (e) the angular momentum increases due to the numerical 
error for Fi = 1.325 but decreases for Fi = 1.3. The reason for (a) is clear because for smaller values of Fi, the 
central density increases by a large factor, and as a result, the grid resolution becomes worse even in the identical grid 
spacing. The reason for (b) is also clear because with stiffer equations of state, the central density of neutron stars 
is smaller. The result (c) comes from the fact that with smaller values of Fi, shocks are stronger, and as a result, 
the amount of mass that is ejected outside the core is larger. However, the reason for (d) and (e) is not clear at all. 
These may be consequences of our choice of the finite differencing scheme for the hydrodynamic equations. Indeed, a 
similar tendency is found in the results of an approximate general relativistic simulation [14]. 

It is interesting to note that even for (Clb), in which the maximum allowed baryon rest-mass of neutron stars is 
very small (~ IMq), a low-mass neutron star of M^core ~ O.55M0 is formed after the collapse. This is because the 
shocks explode the infalling matter sufficiently. Thus, even in the case that the stellar core mass before collapse is 
much larger than the maximum allowed mass of the neutron star for a given equation of state, a neutron star instead 
of a black hole could be formed at least temporarily, although such a neutron star could easily collapse to a black hole 
by subsequent accretion of matter or by a fall-back. 

In Fig. 14, the specific angular momentum spectra at selected time steps are shown for (Cld) and (C2a) as examples. 
In both cases, we take N — 600. The figures demonstrate that the spectral shape is well conserved throughout the 
simulations irrespective of the initial angular velocity profile. Therefore, we conclude that angular momentum transfer 
due to numerical dissipation is sufficiently small in numerical computations. 



D. An example of stellar core collapse 

In this subsection, we present numerical results for a simulation of rotating stellar core collapse that is started 
from a realistic initial condition with the central density ~ 10^*^ g/cm^ and with F = 4/3. We give a rigidly rotating 
equilibrium star that is close to the mass shedding limit as the initial condition. Several quantities of this rotating 
star arc listed in Table V. Since it is rapidly rotating and its compactness is sufficiently small, this equilibrium star 
is dynamically stable against gravitational collapse even for the polytropic equation of state with F = 4/3 [43]. Thus, 
the collapse is triggered by the slight decrease of the adiabatic constant from 4/3 to Fi. In the present simulation, 
we set Fi = 1.325 and F2 = 2. 
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(a) (b) 
FIG. 12. The same as Fig. 11 but for collapse of a star (C2a). In both figures, the solid, dotted, and dashed curves denote 
the results with N = 600, 400, and 300. With these grid numbers, Aa; =0.385, 0.577, and 0.770 km. For t > 7 msec, shocks 
reach outer boundaries of the computational domain, and matter starts escaping. 



Pc (g/cm'') 


M^Mq) 


M(M0) 


R (km) 


\T/W\ 


J/M'-' 


ac 


1.65 X 10'" 


1.491 


1.491 


1910 


8.89 X lO-" 


1.136 


0.993 



TABLE V. Central density, baryon rest-mass, ADM mass, equatorial radius, ratio of the rotational kinetic energy to potential 



energy, non-dimensional angular momentum parameter, and central value of the lapse function of a rotating star chosen as an 
initial condition for a stellar core collapse simulation in Sec. IV D. 



Since the characteristic length scale changes by a factor of ~ 100 during the collapse, wc performed the simulation 
changing the grid size and grid number as done in [22]. The grid size and computational domain were changed 
monitoring the value of the lapse function at the center (ac), which approximately indicates the compactness of the 
collapsing star. Whenever wc carried out rcgridding, we made the grid spacing half and used cubic interpolation 
[32] for assigning the values of variables on the finer grids. The simulation was started with N = 500, by which the 
equatorial radius is covered by 480 grid points initially. At t = 0, ac ~ 0.993. We carried out the first regridding when 
ttc was 0.975, at which the mean radius of the collapsing star became ^ 1/4 of the initial one. In this rcgridding. wc 
chose A'^ = 900 and made the grid spacing half. The next rcgridding was carried out when ac = 0.95 and 0.90, and 
we chose A'' = 1500 and 2100, respectively. After ac reached 0.90, we fixed N and grid spacing. L and Ax in the final 
stage are about 1050 km and 0.5 km, respectively. For this simulation, the computational time was about 100 CPU 
hours for « 60000 time steps using eight processors of the FACOM VPP 5000 machine. 

Since the computational region was reduced whenever we carried out the regridding, a small amount of mass that is 
outside the new computational domain was discarded. However, the magnitude of the violation of mass and angular 
momentum conservation is less than ~ 0.5% and 2%, respectively [see Fig. 16(b)]. This implies that the total amount 
of the discarded mass is comparable to that of the numerical error associated with the finite differencing, which does 
not much affect the evolution of the system, as indicated in Sec. IV C. 

In Fig. 15, we display the density contour curves and velocity fields at selected time steps around which shocks are 
formed. The time of the shock formation is ^ 71.7 msec, which is in good agreement with that for model A1B3G1 in 
[14] with the correction factor which is associated with the dynamical time scale as (/3c,init/10^° g/cm^)"-'^/^. (Note 
that in [14], the central density of the initial condition is 10^" g/cm"^, while here Pc,imt ~ 1-65 x 10^° g/cm^.) This 
coincidence suggests that the approximate general relativistic approach adopted in [14] is indeed suited for study of 
axisymmetric stellar core collapse to neutron stars in general relativity. 

After the shock formation, the shock fronts of prolate shape spread outward. The prolateness is produced by the 
fact that the shocks are stronger for the z direction due to the absence of centrifugal force [11,12,14]. In Fig. 16(a), 
we display the time evolution of the central density and lapse function at the center. Global features are qualitatively 
the same as those for the simulation of model (Cla) presented in Sec. IV C. As in that case, the central density (lapse 
function) monotonically increases (decreases) until it exceeds Pnuc- When it becomes ~ 3.5/9„uc, the collapse is halted 



20 





02468 0246 
T (msec) T (msec) 



(a) (b) 
FIG. 13. (a) Time evolution of central density, lapse function at origin, and fraction of baryon rest-mass for p > pnuc, 
and (b) time evolution of violation of rest-mass conservation, angular momentum conservation, and averaged violation of 
the Hamiltonian constraint for (Cla) (solid curves), (Clb) (dotted curves), (Clc) (dashed curves), and (Old) (dotted-dashed 
curves). We adopt N = 600, by which the equatorial radius is initially covered by N grid points for all the models. 
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FIG. 14. Specific angular momentum spectra at selected time steps (a) for (Cld) and (b) for (C2a). 



and the shocks start propagating outward, while the core gradually settles down toward a quasistationary state of 

Pc ^ 2pnuc- 

Although these qualitative features are the same as those for (Cla), there arc a few quantitative differences between 
the results of two models. First, the maximum density and central density of the formed neutron star found here are 
slightly smaller than those for (Cla). This is likely due to the fact that the effect of the centrifugal force plays a more 
important role than for (Cla). Second, the core does not relax to a static state soon, but oscillates approximately in a 
quasiperiodic manner for several periods. This oscillation is not conspicuous for model (Cla). These facts imply that 
to obtain quantitative outputs of stellar core collapse, we should start the simulations from an initial condition of a 
realistic density profile, although qualitative global features of the collapse can be found even using a more compact 
initial condition. We note that the quasiperiodic oscillation found here is also observed in [14]. As reported in [12,14], 
quasiperiodic gravitational waves are likely emitted associated with this oscillation. However, we have not tried to 
compute gravitational waves in the present work. As expected from the results in [14], the amplitude of gravitational 
waves is not very large, so that it would not be technically easy to extract them from the metric in which gauge modes 
and numerical noises are included. Developing a method for the wave extraction of a weak signal will be one of the 
challenging problems in the future. 

In the last figure of Fig. 16(b), the time evolution of baryon rest-mass of the core (of density larger than pnuc) is 
shown. As in the time evolution of pc, it reaches the maximum at which pc reaches the maximum value, and then 
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FIG. 15. Density contour curves of p and velocity fields of at selected time steps around which shocks are formed. The 
contour curves are drawn for p/pnuc = 10~°'*^ , for j = 0, 1, 2, • • • , 20. 



settles down toward a constant ~ O.55M0. Thus, the temporary product after the collapse is a low mass neutron star 
in this simulation. 
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FIG. 16. (a) Time evolution of the central density and lapse function at the center, (b) Time evolution of the violation of 
rest-mass conservation, angular momentum conservation, and the time evolution of rest-mass of the core with p > pnuc- 



V. SUMMARY 



We have presented numerical results obtained by an axisymmetric general relativistic implementation, and demon- 
strated that with this new implementation, it is feasible to carry out long-term simulations for spherical and rapidly 
rotating neutron stars, and rotating stellar collapses to a neutron star and a Kerr black hole. It is shown that the 
simulations for stable neutron stars can be continued for > 10 dynamical time scales until the crash of computation 
with TV ~ 200, even if the neutron stars are of ^ 95% of the maximum allowed mass. The duration of the compu- 
tation is long enough to obtain oscillation modes of neutron stars. The simulations are also feasible for collapse of 
rotating neutron stars to Kerr black holes. We have illustrated that with our implementation, the mass of a Kerr 
black hole formed after the collapse can be computed accurately. We have also demonstrated that it is feasible to 
perform the simulations of rotating stellar core collapse to a neutron star, adopting parametric equations of state that 
mimic realistic equations of state. This illustrates that the new implementation works well also for realistic equations 
of state that have not been adopted so far in fully general relativistic simulations. In conclusion, the axisymmetric 
numerical implementation presented here will be used for a wide variety of astrophysical simulations such as rotating 
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core collapse of a massive star to a neutron star or a black hole, and accretion-induced collapse of a neutron star to a 
black hole. As the next step, we plan to perforin simulations for rotating core collapse to neutron stars, and compare 
the results with those in an approximate relativistic approach [14] . We also consider that collapse of a massive stellar 
core to a black hole is one of the most interesting topics. 

Because of the assumption of axial symmetry, we were able to carry out a wide variety of tests and calibrations for 
our new hydrodynamic implementation with low computational costs. In previous works, e.g., [18,28], tests of their 
hydrodynamic implementations picking up single stars have been done in three spatial dimensions. In those cases, N 
could be at most 100, since the computational costs were very high for N > 100. It is pragmatically very difficult 
to investigate the accuracy and convergence in a well-resolved simulation with N ^ several hundreds under normal 
circumstances in which we can use at most ~ 1000 CPU hours per year. In the axisymmetric case, the simulation 
with A'' ~ several hundreds is not very expensive, and thus for detailed tests of new hydrodynamic implementations 
in general relativity, axisymmetric simulation has great advantages. We expect that for the testing of new gauge 
conditions and wave extraction techniques, it would also play an important role. 

Finally, wc note the following point. Although we focus only on the axisymmetric simulation in this paper, the 
present hydrodynamic implementation can be used for three-dimensional simulations with a slight modification, since 
the transport terms in the hydrodynamic equations are of the same form. Actually, we have already implemented it 
and checked that it works. We expect that with the same grid resolution that we adopted in this paper, the same 
results will be obtained (although it takes a much longer time to carry out the long-term simulations). So far, we have 
performed simulations for a merger of binary neutron stars only using the F-law equation of state [19,20]. However, 
with the new hydrodynamic implementation reported in this paper, wc will be able to adopt a variety of equations 
of state. We plan to perform the simulation for a merger of binary neutron stars adopting more realistic equations of 
state in the future. 
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APPENDIX A: TREATMENT FOR TRANSPORT TERMS IN THE HYDRODYNAMIC EQUATIONS 

Equations (26)-(29) are of the forms 

dtQa + dAF^ = Sa, (Al) 

where Qa and for a = 1-5 are defined as 

Qa = {P*, Jx, Jy. Jz. E,), (A2) 

Fa - [P*v'^, Jxv^ + P»V^Si, Jyv^, J,v^ + Pa^Si, E,v^ + P^{v^ + /3^)]. (A3) 

Here, Ji = p»Ui, E^, = p*e, and Sa in Eq. (Al) denote the right-hand sides of Eqs. (26)-(29). We note that 7 here is 
the determinant of the three- metric in the Cartesian coordinates. 

In numerical computation, wc evaluate the transport terms using the approximate Riemann solver, which relies on 
the characteristic decomposition of the equations (see e.g. [25,27] and references therein). To adopt this method, we 
first need to compute the Jacobian matrix and then to carry out the spectral decomposition of it. 

The Jacobian matrix for the A-ih. direction, M^, is defined by 

{A = xorz). (A4) 

Using this, Eq. (Al) may be expressed in the form 

5 

dtQa + Y. E MabdAQb = Sa. (A5) 
6=1 A=x^z 
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Thus, has information on the characteristic speed of the fluid. 
Following Font et al. [23,24], we calculate the Jacobian matrix from 



^ dqc dQb 
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(A6) 
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where 
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1 dP 
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-fij — ^ij + 

hi = l + e + x 

h2 = 1 + K. 
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B^f^ is obtained by appropriate exchanges of subscripts among x, y, z of B^f^. 

The eigenvalues of the matrix Af^, A"*, correspond to the characteristic speeds of the fluid in the A-th direction, 
and are derived from the equation 



det{B^b - ^^Cab) = 0. 



(A17) 



The solutions are [24,25] 
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= A^, (triple), 



(A18) 



where 



1 

± a2 - VkV'cl 



and 



v^a'il-cD-P^lia'-VkV") 

± acs^{a^ - VkV^){-i^^{a^ - VkV^cfj - (1 - cl)V^V^] 
(no summation for A) 



(A19) 



(A20) 

y'^ = = + /J'^. (A21) 
Using the eigenvalues, the spectrum decomposition for can be done in a straightforward manner as 

M^, = Y,RiXc{R'')ea^ (A22) 

b,c 

where is the diagonal matrix composed of A"^ in the following order: [\^,v^, ,X^]. For convenience of the 

calculation of R^^, we define a matrix T^, which satisfies the relation as 

5 

Rtb = Y.^acTl (A23) 

c=l 

Since i?^^ is calculated from the right eigenvectors of M^, is composed of vectors t^^^ that satisfy the equation 



^(Sf, - X^Ca,){ttf''^ =0 for / = 1 ~ 5. 



(A24) 



6=1 



Then, = [(t^^)(i), (if (tA)(3)^ (^A)(4)^ (^A)(5)]^ and hence we obtain 
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where 



if^'=(A) = 
iJ^'=(A) = 



-c^(«^ - A){q^7^*^ - y''(/3^ + A)} 
-cl{v' - X){a^Y-^ - V'^iP'^ + A)} 



(A27) 
(A28) 



pw'^{v^ — A)^ 

After the above reconstruction of the fiuid equation, the numerical fiuxes in the upwind scheme are computed from 
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(A29) 



where we omit the subscripts for Rab and Kab- Qc and denote Qc at the left and right sides of the corresponding 
interfaces, and are evaluated using the third-order spatial interpolation. At the interface between the i-th and (i-l-l)-th 
cells, we define them according to 

Qi = Qi + ^ + ^, (A30) 
D o 

where Aj = Qi+i — Qi- To suppress the oscillation near shock discontinuities, we modify the interpolation using the 
following min-mod limiter as [46] 

where 

4 = Ai+i/Ai, (A34) 

r- = A,_i/A„ (A35) 

$(r) = minmod(l, br) (1 < 6 < 4 for TVD condition). (A36) 

For the simulations presented in this paper, we choose 6 = 2, since for fe « 1, the dissipation is so large that the 
envelope of neutron stars spreads outward too quickly, while for 6 = 4, the oscillation around shock discontinuities is 
too serious. 

The values for components of matrices Rab and Aab at grid interfaces are computed using the Roe-type average 
such as [45,46] 



<}i+l/2 = y===-—== , (A37) 

where we carry out the average for variables Ui, k, x, and h (i.e., qi is one of these variables). Other variables are 
computed from them. It should be noted that in the relativistic case, the average is not uniquely specified in contrast 
with the Newtonian case [45]. However, numerical results of test computations (see below) seem to indicate that this 
averaging is appropriate. 

To confirm that our hydrodynamic implementation can capture shocks accurately, we carried out the simulations 
for Riemann shock tube problems and wall shock problems in the 1-1-1 special relativistic spacetime with {t,x) as the 
coordinates. In this test, we adopt the F-law equation of state. In both tests, we take TV = 400 with Aa; = 1/A''. 

In the Riemann shock-tube problem, we choose T = 5/3. The parameters of the initial condition are chosen as 
p = 10 and P = 13.3 for a; < and p = I and P = 10~^ for a; > following previous papers [27]. In the wall shock 
problem, we set the parameters as t;^ = 0.9, p= 1, and P = 10~^ with T = 4/3. 

In Figs. 17 and 18, we compare numerical results of the Riemann shock tube problem and of the wall shock problem 
with analytical solution (solid curves) for the choice 6=1 and 2. As indicated in these figures, numerical results agree 
well with analytic solutions, in particular for 6 = 2. 
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